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A method is presented for the unbiased numerical computation of two-particle response functions 
of correlated electron materials via a solution of the dynamical mean-field equations in the presence 
of a perturbing field. The power of the method is demonstrated via a computation of the Raman 
Big and E>2 g scattering intensities of the two dimensional Hubbard model, in parameter regimes 
believed to be relevant to high-temperature superconductivity. The theory reproduces the 'two- 
magnon' peak characteristic of the Raman intensity of the insulating parent compounds of high-T c 
copper oxide superconductors and shows how it evolves to a quasiparticle response as carriers are 
added. The method can be applied in any situation where a solution of the equilibrium dynamical 
mean-field equations is feasible. 
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The development of dynamical mean-field theory, first 
in its single-site version [1 and then in cluster exten- 
sions 2 , along with its interface to band theory [3J 3] , 
has transformed our understanding of correlated electron 
physics. In the particular case of the two-dimensional 
Hubbard model, believed [5] to represent the physics of 
high temperature superconductors, the method has pro- 
vided new insights into the correlation-driven ('Mott') 
metal-insulator transition (HJ |7j, the pseudogap regime 
that separates the Mott insulator from the Fermi liquid 
metal [7HT0] and the existence [1TT - [T3] and properties [T4l - 
117] of a d-wave superconducting state. However, dynam- 
ical mean-field theory is based on an approximation to 
the one-electron Green function G(k,u>), measurable in 
angle-resolved photoemission experiments |18| while the- 
oretical analysis of wide classes of experiments including 
optical conductivity, Raman spectroscopy and inelastic 
neutron scattering requires a vertex function whose com- 
putation has proven very challenging. While the charge 
vertex can be obtained analytically in simplified situa- 
tions such as the Falicov-Kimball model 19J, computa- 
tions of the vertex for interacting electron models have 
not, in practice, been carried out in full generality. Kuncs 
|20j obtained the zero frequency charge and spin vertices 
corresponding to single-site dynamical mean-field theory 
of the Hubbard model, and Yang et al. [21] obtained 
the spin, charge and superconducting vertices for larger 
clusters. However, calculations [551 123] of the full dy- 
namic (uj 0) response have been based on the single 
site approximation and have employed truncations of the 
general frequency dependence. 

In this paper we present a new method for determining 
the two-particle response in cluster dynamical mean-field 
theory and demonstrate its effectiveness via a computa- 
tion of the doping dependence of the Raman scattering 
amplitude of the two-dimensional Hubbard model from 
the Mott insulating to the Fermi liquid regime. Raman 



spectroscopy has been of fundamental importance to high 
temperature superconductivity |24j but the theoretical 
description in terms of an underlying Hubbard model in- 
volves vertex corrections in an essential way and has not 
been systematically studied. 

To introduce our method we recall salient features of 
the theory of linear response |25] , defined quantum me- 
chanically as the leading-order difference of the expec- 
tation value of an operator R in the presence and ab- 
sence of a probe field P. This is given by (R) p = 

Tr [rG p ] - Tr [RG e9 j = xupP + 0{P 2 ) with 
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Here G cq = (Gq 1 - S^)" 1 is the equilibrium (P = 0) 
Green function, related in the standard way to a bare 
Green function Go and a self energy X. We have omit- 
ted a possible term arising from explicit dependence of 
R on P; this gives rise to the 'diamagnetic' term in the 
optical conductivity but is not otherwise relevant. The 
first term in Eq. [T] gives the 'bubble term', computable 
from knowledge of the one-electron Green function and 
the bare vertex SGq 1 /SP; the second term, arising from 
changes in the many-body physics due to the perturba- 
tion, is the vertex correction of interest here. 

For wide classes of strongly correlated materials, 
neither perturbative diagrammatic expansions about a 
mean-field solution nor partial (e.g RPA or GW) resum- 
mations suffice; a fully nonperturbative treatment is re- 
quired. For the one electron Green function, cluster dy- 
namical mean- field theory [5] provides such a treatment. 
In this theory the electron self energy S, a matrix in 
the full single-particle Hilbert space of the problem, is 
approximated in terms of a much smaller number N c of 
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Different choices for (f> a p give different versions of dynam- 
ical mean-field theory while the T, a/3 are the self energies 
of a quantum impurity model specified by the interac- 
tions of the original model and by mean-field functions 
(Gq 1 ) 01 ^ determined by the self consistency condition 
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Here the trace and inner inversion are over the single 
particle Hilbcrt space of the full problem and the outer 
inversion is in the impurity model space. All quanti- 
ties may depend on the perturbation P and we compute 
SH/SP directly by solving to first order in P. 

Linearization of Eq. [3] yields the first order correc- 
tion (Gq 1 )^^ in terms of SGq 1 /SP (known a-priori) and 
8Y, afs /8P (to be computed). Fr om this and the im- 
purity model four-point function T, written here for a 
monochromatic perturbation of frequency Q as 

r{?#Vo/) = ( c jc + f7)4H v ( w ')cL( w ' + o)) (4) 

we obtain the first order correction G 1 to the impurity 
model Green function as 

-fe.oG c Q >)£ ^pMiGo^pfa'Ui)' ( 5 ) 

wi ,ct' (3' 

From this and the linearized impurity-model Dyson equa- 
tion follows S"£ a P /5P. The resulting linear equation is 
solved for <5£ Q ' 3 /SP, from which we obtain the desired 
vertex correction SH/SP via Eq. [5] 

Previous dynamical mean-field literature introduced 
[HUH, and used [20J [211231 EH], a different approach, in- 
verting the impurity model Bethe-Saltpeter equation to 
obtain the two-particle irreducible impurity-model vertex 
in terms of T and then using this in the lattice Bethe- 
Saltpeter equation to compute the physical vertex. Our 
procedure replaces the numerical inversion of the Bethe- 
Saltpeter equation (which requires considerable care [20] 
to avoid numerical instabilities [221 [271 EE] ) by the solu- 
tion of the linearized DMFT equation (which we have 
found not to be problematic), and avoids the second 
Bethe-Saltpeter equation by constructing the fully re- 
ducible lattice vertex directly. 

The key issue in the implementation is the measure- 
ment and storage of T, which is needed in each sector 
for a wide range of u>, u>\ at every relevant Q. It is 
necessary to compute T in a strip \oj — u>i\ < A\ and 
\ui + oji + fi| < |f2| + A2. Ai 2 are increased empirically 



until the final \ ceases to change on the ~ 1% level and 
we typically need Ai ~ A2/2 of order 3-4 times a relevant 
frequency scale (interaction strength or bandwidth) . T is 
measured on the imaginary time axis and the Fourier 
transform to frequencies is a significant portion of the 
computational burden. Computations are substantially 
accelerated by a formal rearrangement which allows the 
needed information to be obtained from a small set of one 
dimensional Fourier transforms rather than a two dimen- 
sional one, and by recent extensions of the fast Fourier 
transform method to non-uniform grids |29) . 
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FIG. 1: Raman B\ g scattering intensity in the 4-site (solid 
line, red online) and 8-site (dashed line, green on-line) DCA 
approximations for parameters U = It, t! = — 0.15t, f3 = 10/t 
and n — 1. Upper right inset: interaction strength depen- 
dence of Big Raman intensity at density n = 1 calculated 
for the four site cluster. Lower two insets: momentum space 
partitioning for the four- and eight-site cluster geometries con- 
sidered in this paper. 



To demonstrate the power of the new approach we 
compute the non-resonant Raman scattering intensity of 
the two dimensional Hubbard model defined by 



#Hub = H + U ^ "iT n 4- 
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with 
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For definiteness we take Hq = 
£fe = — 2t (cosk x + cos k y ) 
= —0.154 and U — It. To fix the energy scale we use 
t = 0.35 eV, a value generally accepted for high-T c su- 
perconductors [50], 

We consider two scattering geometries: P>i g , where the 
electric fields of the incident and outgoing photons are di- 
rected along the Brillouin zone axes {k x or k v = 0) high- 
lighting the antinodal region, and B2 g , where the electric 
field vectors are directed along Brillouin zone diagonals 
(k x — ±k y ) highlighting the nodal region. The perturb- 
ing terms corresponding to the B\ g and B>i g scattering 
channels are (in the non-resonant approximation) 
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We use the dynamical cluster approximation (DCA) 
[2] in which the Brillouin zone of momentum space is 
partitioned into N c equal area tiles labeled by central 
momentum K and study in particular N c = 4 and 8 (see 
insets to Fig. [I]). Both lattice and impurity E in Eq. [2] 
are diagonal. The lattice quantities depend on a contin- 
uous momentum k, the index a = /3 represents cluster 
momentum K and § k K = 1 if k is in tile K and zero 
otherwise so the lattice trace is just a momentum inte- 
gral over the tile. We solve the model on the imaginary 
axis in the paramagnetic phase using the numerically ex- 
act continuous-time auxiliary field (CT-AUX) impurity 
solver [31j with submatrix updates [3 2) and analytically 
continue the final XRp(^m) using the maximum entropy 
technique 33J taking into account covariance matrices es- 
timated by a jackknife procedure. 
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FIG. 2: Raman B\ g scattering function of isolated 4-site clus- 
ter calculated at U = St and ft — 20 jt by exact diagonaliza- 
tion (solid lines) and Quantum Monte Carlo (dashed lines) 
for Matsubara (main panel) and real frequencies (inset). 

To test the formalism we observe that if the momentum 
integral is replaced by an evaluation at the central mo- 
mentum K then our procedure is reduced to the solution 
of an isolated cluster, which can also be solved by exact 
diagonalization (ED). The main panel of Fig. [2] compares 
the Matsubara axis Raman B\ g scattering intensity ob- 
tained by applying our procedure to the isolated four site 
cluster with that obtained from a direct diagonalization 
of the same isolated cluster. The results are seen to be 
identical up to Monte Carlo errors which are smaller than 
the symbol size. The inset compares the real axis spec- 
trum obtained by analytical calculation to that obtained 



directly from the exact solution. Apart from a broaden- 
ing, the two procedures give the same result; in particular 
the areas of the peaks are the same in the two methods. 

Fig. [T] shows the calculated Raman B lg intensity at 
carrier density n = 1, i.e. in the paramagnetic Mott in- 
sulating phase. The spectra exhibit a two peak structure. 
The peak at higher energies corresponds to quasiparticle 
excitations across the Mott gap ~ uj — 3t for the 8-site 
(integrated area 2.5) and ~ ui — 5t for the 4-site ap- 
proximation (integrated area 3.1). We identify the lower 
energy peak (area 1.0 for 4-site and 0.2 for 8-site) as 
arising from the creation of a pair of spin flip excitations 
because at n = 1 and large U these are the only excita- 
tions available in this energy range. This identification 
is corroborated by computations (inset) of the interac- 
tion strength dependence in the 4-site cluster (present- 
day computational limits preclude study of U > 7t for 
the 8-site cluster): as U is increased the upper feature 
shifts up in energy and the lower feature shifts down, as 
expected for a peak scaling with J ~ t 2 /U. 

While the 4 and 8 site calculations are qualitatively 
similar, they are quantitatively different. In the DCA ap- 
proximation the 4-site cluster is known [7j to have prop- 
erties different from all of the other clusters, among other 
things strongly overestimating singlet formation and in- 
sulating behavior. In the rest of this paper we focus on 
the 8-site cluster, believed [7] to be more representative 
of the physics of the model. 
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FIG. 3: Main panel: Raman B\ g scattering intensity in 8-site 
cluster for parameters U = 7t, t' — — O.lbt and /3 = 20 /t at 
dopings indicated. Inset: separation into bubble and total 
contribution for x — 0.065. 

Fig. [3] presents the evolution of the B lg spectra. When 
the insulator is doped away from half-filling, the peak at 
u) ~ t broadens and shifts to lower frequency, while at 
lowest frequencies a component x( w ) ~ w appears and 
increases rapidly with doping. The inset decomposes the 
x = 0.065 spectrum into bubble and vertex correction 
(visible as the difference between the two curves). The 
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bubble diagram accounts for the entire quasiparticle part 
while the vertex correction, which in the insulating state 
gives rise to the two-spin-flip peak, produces the higher 
frequency maximum. We therefore attribute the higher 
frequency peak to the relic of the two spin-flip peak in 
the metallic state and the low frequency x(w) ~ cj feature 
to quasiparticles. At the lowest dopings x — 0.065 and 
x = 0.086 the pseudogap is visible as a change in slope 
from the very low frequency regime (dominated by quasi- 
particles) to an intermediate energy regime where much 
of the scattering comes from the two spin-flip feature. 
As the doping is increased beyond x = 0.1 the pseudogap 
disappears, the two-spin flip and quasiparticle scatterings 
merge, and the vertex correction decreases in importance, 
becoming completely negligible by doping ~ 0.24. These 
data are in excellent semiquantitative agreement with the 
Big spectra reported in Fig. 9 of Ref. 1541 
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FIG. 4: Temperature dependence of Raman B\ g scattering 
intensity in 8-site DCA at hole dopings indicated. 

Fig. [4] shows the temperature dependence of the calcu- 
lated Raman scattering intensity in B\ g channel at dop- 
ing level x = 0.065 (upper panel) and x = 0.16 (lower 
panel). At the lower doping the low-frequency Raman 
Big intensity is suppressed as the temperature decreases, 
whereas at the higher doping the initial slope is seen to 
increase as the temperature decreases again consistent 
with data 34 37 and with exact diagonalization of small 
clusters in the t — J model [55] ■ 

We next turn to the B^g channel, which highlights 
the zone diagonal region of the Brillouin zone. The 
zone diagonal region is not affected by the pseudogap 
and the quasiparticle velocity is high and may therefore 
also be expected to dominate the optical conductivity 
a. This along with a standard relation for metals moti- 
vated Ref. [59"] to identify Irnx.B 2 (w)/w with Rea(ui) [55] , 
The main panel of Fig. [5] presents the Big spectra as 
Imxs 2 (u>)/uj. We see a 'Drude' peak centered at zero 



frequency which grows noticeably and sharpens as dop- 
ing is increased, and a broad higher-frequency continuum 
which is only weakly doping dependent. The inset shows 
that as temperature is varied the 'Drude' peak increases 
in height and decreases in width; there is also a small 
~ 10% (not shown) increase in area with decreasing T. 
These features bear a strong qualitative resemblance to 
the optical conductivity data taken in the high-T c cuprate 
superconductors [4"U1 14"1"] . 
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FIG. 5: Main panel: Big Raman spectra calculated at U = It 
and inverse temperature j3 = 20/t ~ 200 K for hole dopings 
indicated and presented as Im.XB 2g /u. Inset: temperature 
dependence of ImxB 2g / ljJ f° r hole doping x = 0.13. 

In this paper we have developed a new method for 
treating the vertex corrections which are essential for 
the computation of wide classes of experimentally rel- 
evant spectroscopies of interacting electron systems. We 
used the formalism to show that the two dimensional 
Hubbard model accounts for the essential features of the 
doping-dependent Raman spectra observed [34j in high- 
ly copper-oxide superconductors. The formalism intro- 
duced here is straightforwardly generalizable to most 
other response functions, for example to the momentum 
dependent spin response needed for neutron scattering, 
although in the special case of the optical conductivity 
the Ward identity issues discussed in Ref. [32] create com- 
plications which are not yet resolved. The additional 
computational burden of our methods scales as a power 
law in cluster size, whereas the equilibrium dynamical 
mean-field computations themselves are limited by an 
exponential barrier imposed by the fermion sign prob- 
lem and Hilbert space size. Therefore we expect that 
as computers grow more powerful, our methods can be 
applied to essentially any case for which an equilibrium 
DMFT solution is feasible, opening new directions for the 
theoretical understanding of correlated electron materi- 
als. 
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